In this script we are downloading current and future environmental data (as rasters) to use in species distribution models (SDMs) for koalas (Phascolarctos cinereus) in the South-East Queensland (SEQ) region.
Import packages
Code
library(terra)
terra 1.8.50
Code
library(dplyr)
Attaching package: 'dplyr'
The following objects are masked from 'package:terra':
intersect, union
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
Code
library(sf)
Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
Code
library(ggplot2)
Load South East Queensland (SEQ) boundary
We start by defining our study area, which is the South East Queensland (SEQ) region. We will use the Local Government Areas (LGA) shapefile to define the extent of SEQ.
# Load the study area shapefileLGA <-st_read("Data/Environmental_variables/Local_Government_Areas.shp")
Reading layer `Local_Government_Areas' from data source
`/Users/scottforrest/Library/CloudStorage/OneDrive-QueenslandUniversityofTechnology/PhD - Scott Forrest/GIT/ICCB_geospatial_tools_conservation/Session 3/Data/Environmental_variables/Local_Government_Areas.shp'
using driver `ESRI Shapefile'
Simple feature collection with 78 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 137.9946 ymin: -29.17927 xmax: 153.5519 ymax: -9.087991
Geodetic CRS: GDA94
Code
# Check the coordinate reference system (CRS)st_crs(LGA)
Coordinate Reference System:
User input: GDA94
wkt:
GEOGCRS["GDA94",
DATUM["Geocentric Datum of Australia 1994",
ELLIPSOID["GRS 1980",6378137,298.257222101,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore."],
BBOX[-60.55,93.41,-8.47,173.34]],
ID["EPSG",4283]]
Code
# Convert to WGS84LGA <- LGA %>%st_transform(3112)# Select local govt. areas for South East QueenslandLGA_SEQ <- LGA %>%filter(lga %in%c("Brisbane City", "Moreton Bay City", "Logan City", "Ipswich City", "Redland City", "Scenic Rim Regional", "Somerset Regional", "Lockyer Valley Regional", "Gold Coast City", "Sunshine Coast Regional", "Toowoomba Regional", "Noosa Shire"))ggplot() +geom_sf(data = LGA, color ="black") +geom_sf(data = LGA_SEQ, fill ="purple3", alpha =0.5, color ="black", size =0.2) +theme_minimal() +theme(legend.position ="none") +labs(title ="Local Government Areas Queensland (SEQ in purple)")
Code
ggplot() +geom_sf(data = LGA_SEQ, fill ="purple3", alpha =0.5, color ="black", size =0.2) +theme_minimal() +theme(legend.position ="none") +labs(title ="Local Government Areas South East Queensland (SEQ)")
Merge into a single polygon
Code
# Merge the SEQ LGAs into one polygonSEQ_extent <-st_union(LGA_SEQ)ggplot() +geom_sf(data = SEQ_extent, fill ="purple3", alpha =0.5, color ="black", size =0.2) +theme_minimal() +theme(legend.position ="none") +ggtitle("South-East Queensland Spatial Extent") +theme_bw()
Save SEQ extent for other scripts
Code
# Convert our SEQ extent to a SpatExtent object by converting to a SpatVectorSEQ_extent.vect <- terra::vect(SEQ_extent)writeVector(SEQ_extent.vect, "Data/Environmental_variables/SEQ_extent.shp", overwrite = T)
Load current environmental data
Layers were made available to us by the EcoCommons team and were created by Toombs and Ma (2025):
Toombs, N., and Ma S., 2025, A High-Resolution Dataset of 19 Bioclimatic Indices over Australia, Climate Projections and Services – Queensland Treasury, Brisbane, Queensland. [https://longpaddock.qld.gov.au/qld-future-climate/data-info/tern/]
Code
files <-list.files("Data/Environmental_variables/Current_climate_QLD", pattern =".tif$", full.names =TRUE)# Load all bioclim rasterscurrent_bioclim <-lapply(files, terra::rast) # Make into one raster stackcurrent_bioclim <-rast(current_bioclim)plot(current_bioclim)
current_bioclim <- terra::mask(current_bioclim, SEQ_extent.vect)# You can see that this has masked the area but the extent is still the sameplot(current_bioclim[[1]])
# Save the current environmental covariateswriteRaster(current_bioclim, filename ="Data/Environmental_variables/current_bioclim.tif",overwrite = T)
Load future environmental data
Here we load outputs from a moderate-high emissions shared socio-economic path scenario (SSP 3.70) for the year 2090 (2080 - 2099).
Code
files <-list.files("Data/Environmental_variables/Future_climate_SSP370_2090", pattern =".tif$", full.names =TRUE)# Load all bioclim rastersfuture_bioclim <-lapply(files, terra::rast) # Make into one raster stackfuture_bioclim <-rast(future_bioclim)plot(future_bioclim)
future_bioclim <- terra::subst(from =0, to =NA, future_bioclim) # Set all values of 0 to NAfuture_bioclim <- terra::mask(future_bioclim, SEQ_extent.vect)# You can see that this has masked the area but the extent is still the sameplot(future_bioclim[[1]])
# Save the future environmental covariateswriteRaster(future_bioclim, filename ="Data/Environmental_variables/future_bioclim.2090.SSP370.tif",overwrite = T)
Load future environmental data 2
Here we load outputs from a low emissions shared socio-economic path scenario (SSP 1.26) for the year 2090 (2080 - 2099).
Code
files <-list.files("Data/Environmental_variables/Future_climate_SSP126_2090", pattern =".tif$", full.names =TRUE)# Load all bioclim rastersfuture_bioclim <-lapply(files, terra::rast) # Make into one raster stackfuture_bioclim <-rast(future_bioclim)plot(future_bioclim)
future_bioclim <- terra::subst(from =0, to =NA, future_bioclim) # Set all values of 0 to NAfuture_bioclim <- terra::mask(future_bioclim, SEQ_extent.vect)# You can see that this has masked the area but the extent is still the sameplot(future_bioclim[[1]])
# Save the future environmental covariateswriteRaster(future_bioclim, filename ="Data/Environmental_variables/future_bioclim.2090.SSP126.tif",overwrite = T)
---title: "ICCB Environmental data download"author: "Scott Forrest and Charlotte Patterson"date: "`r Sys.Date()`"execute: cache: false# bibliography: references.bibtoc: truenumber-sections: falseformat: html: self-contained: true code-fold: show code-tools: true df-print: paged code-line-numbers: true code-overflow: scroll fig-format: png fig-dpi: 300 pdf: geometry: - top=30mm - left=30mmeditor: sourceabstract: | In this script we are downloading current and future environmental data (as rasters) to use in species distribution models (SDMs) for koalas (Phascolarctos cinereus) in the South-East Queensland (SEQ) region. ---## Import packages```{r}library(terra)library(dplyr)library(sf)library(ggplot2)```## Load South East Queensland (SEQ) boundaryWe start by defining our study area, which is the South East Queensland (SEQ) region. We will use the Local Government Areas (LGA) shapefile to define the extent of SEQ.https://qldspatial.information.qld.gov.au/catalogue/custom/detail.page?fid={3F3DBD69-647B-4833-B0A5-CC43D5E70699}```{r}# Load the study area shapefileLGA <-st_read("Data/Environmental_variables/Local_Government_Areas.shp")# Check the coordinate reference system (CRS)st_crs(LGA)# Convert to WGS84LGA <- LGA %>%st_transform(3112)# Select local govt. areas for South East QueenslandLGA_SEQ <- LGA %>%filter(lga %in%c("Brisbane City", "Moreton Bay City", "Logan City", "Ipswich City", "Redland City", "Scenic Rim Regional", "Somerset Regional", "Lockyer Valley Regional", "Gold Coast City", "Sunshine Coast Regional", "Toowoomba Regional", "Noosa Shire"))ggplot() +geom_sf(data = LGA, color ="black") +geom_sf(data = LGA_SEQ, fill ="purple3", alpha =0.5, color ="black", size =0.2) +theme_minimal() +theme(legend.position ="none") +labs(title ="Local Government Areas Queensland (SEQ in purple)")ggplot() +geom_sf(data = LGA_SEQ, fill ="purple3", alpha =0.5, color ="black", size =0.2) +theme_minimal() +theme(legend.position ="none") +labs(title ="Local Government Areas South East Queensland (SEQ)")```## Merge into a single polygon```{r}# Merge the SEQ LGAs into one polygonSEQ_extent <-st_union(LGA_SEQ)ggplot() +geom_sf(data = SEQ_extent, fill ="purple3", alpha =0.5, color ="black", size =0.2) +theme_minimal() +theme(legend.position ="none") +ggtitle("South-East Queensland Spatial Extent") +theme_bw() ```## Save SEQ extent for other scripts```{r}# Convert our SEQ extent to a SpatExtent object by converting to a SpatVectorSEQ_extent.vect <- terra::vect(SEQ_extent)writeVector(SEQ_extent.vect, "Data/Environmental_variables/SEQ_extent.shp", overwrite = T)```## Load current environmental dataLayers were made available to us by the EcoCommons team and were created by Toombs and Ma (2025):Toombs, N., and Ma S., 2025, A High-Resolution Dataset of 19 Bioclimatic Indices over Australia, Climate Projections and Services – Queensland Treasury, Brisbane, Queensland. [https://longpaddock.qld.gov.au/qld-future-climate/data-info/tern/]```{r}files <-list.files("Data/Environmental_variables/Current_climate_QLD", pattern =".tif$", full.names =TRUE)# Load all bioclim rasterscurrent_bioclim <-lapply(files, terra::rast) # Make into one raster stackcurrent_bioclim <-rast(current_bioclim)plot(current_bioclim)# Examine the resolutioncurrent_bioclim# Check the CRScrs(current_bioclim)# Update CRS current_bioclim <- terra::project(current_bioclim, "EPSG:3112")# Our resolution is now ~5km by 5kmcurrent_bioclim```## Mask to SEQ extent```{r}current_bioclim <- terra::mask(current_bioclim, SEQ_extent.vect)# You can see that this has masked the area but the extent is still the sameplot(current_bioclim[[1]])current_bioclim <- terra::crop(current_bioclim, SEQ_extent.vect)plot(current_bioclim)# Save the current environmental covariateswriteRaster(current_bioclim, filename ="Data/Environmental_variables/current_bioclim.tif",overwrite = T)```## Load future environmental dataHere we load outputs from a moderate-high emissions shared socio-economic path scenario (SSP 3.70) for the year 2090 (2080 - 2099).```{r}files <-list.files("Data/Environmental_variables/Future_climate_SSP370_2090", pattern =".tif$", full.names =TRUE)# Load all bioclim rastersfuture_bioclim <-lapply(files, terra::rast) # Make into one raster stackfuture_bioclim <-rast(future_bioclim)plot(future_bioclim)# Examine the resolutionfuture_bioclim# Check the CRScrs(future_bioclim)# Update CRS future_bioclim <- terra::project(future_bioclim, "EPSG:3112")# Our resolution is now ~5km by 5kmfuture_bioclim```## Mask to SEQ extent```{r}future_bioclim <- terra::subst(from =0, to =NA, future_bioclim) # Set all values of 0 to NAfuture_bioclim <- terra::mask(future_bioclim, SEQ_extent.vect)# You can see that this has masked the area but the extent is still the sameplot(future_bioclim[[1]])future_bioclim <- terra::crop(future_bioclim, SEQ_extent.vect)plot(future_bioclim[[1]])# Save the future environmental covariateswriteRaster(future_bioclim, filename ="Data/Environmental_variables/future_bioclim.2090.SSP370.tif",overwrite = T)```## Load future environmental data 2Here we load outputs from a low emissions shared socio-economic path scenario (SSP 1.26) for the year 2090 (2080 - 2099).```{r}files <-list.files("Data/Environmental_variables/Future_climate_SSP126_2090", pattern =".tif$", full.names =TRUE)# Load all bioclim rastersfuture_bioclim <-lapply(files, terra::rast) # Make into one raster stackfuture_bioclim <-rast(future_bioclim)plot(future_bioclim)# Examine the resolutionfuture_bioclim# Check the CRScrs(future_bioclim)# Update CRS future_bioclim <- terra::project(future_bioclim, "EPSG:3112")# Our resolution is now ~5km by 5kmfuture_bioclim```## Mask to SEQ extent```{r}future_bioclim <- terra::subst(from =0, to =NA, future_bioclim) # Set all values of 0 to NAfuture_bioclim <- terra::mask(future_bioclim, SEQ_extent.vect)# You can see that this has masked the area but the extent is still the sameplot(future_bioclim[[1]])future_bioclim <- terra::crop(future_bioclim, SEQ_extent.vect)plot(future_bioclim)# Save the future environmental covariateswriteRaster(future_bioclim, filename ="Data/Environmental_variables/future_bioclim.2090.SSP126.tif",overwrite = T)```